# Jay Jang 4/17/2021
# R version 4.0.3 (2020-10-10)
# Platform: x86_64-w64-mingw32/x64 (64-bit)

library(tidyverse)
library(lfe)
library(stargazer)
library(Amelia)
library(xtable)

setwd('C:/Users/battl/OneDrive/JUN/UC MERCED/Independent Research/local election and media freedom')

df <- read_csv('data/clean/df.csv') 


#=======#
# plots #
#=======#

## Figure A1

# local election plot

df %>%
  select(year, v2ellocelc_bi) %>%
  drop_na() %>%
  ggplot(aes(x=factor(year), fill=factor(v2ellocelc_bi)))+
  theme_classic(base_size = 20)+geom_bar(position = "fill")+xlab("year")+ labs(fill = "Local Election")+ylab("Proportion of local election")+
  theme(axis.text.x=element_text(angle=90,hjust=1))
ggsave(filename = "plot/prop_locelec.png",  width = 15, height = 5, dpi = 300, units = "in", device='png') # save last plot

# regional election plot

df %>%
  select(year, v2elsrgel_bi) %>%
  drop_na() %>%
  ggplot(aes(x=factor(year), fill=factor(v2elsrgel_bi)))+
  theme_classic(base_size = 20)+geom_bar(position = "fill")+xlab("year")+ labs(fill = "Regional Election")+ylab("Proportion of Regional election")+
  theme(axis.text.x=element_text(angle=90,hjust=1))
ggsave(filename = "plot/prop_regelec.png",  width = 15, height = 5, dpi = 300, units = "in", device='png') # save last plot


## Figure A2

# density plot of local election - media freedom

df %>%
  select(MSF, v2ellocelc_bi) %>%
  drop_na() %>%
  ggplot(aes(x=MSF, fill=factor(v2ellocelc_bi)))+geom_density(position="identity", alpha=0.5, adjust = 1/2) + theme_bw(base_size = 20) + xlab("Vdem_Media") +
  scale_fill_discrete(name = "Local Election", labels = c("No", "Yes")) +  theme(legend.position = "bottom") + xlab("Media Freedom")
ggsave(filename = "plot/msfs_bi_local.png",  width = 6, height = 5, dpi = 300, units = "in", device='png') # save last plot


# density plot of regional election - media freedom

df %>%
  select(MSF, v2elsrgel_bi) %>%
  drop_na() %>%
  ggplot(aes(x=MSF, fill=factor(v2elsrgel_bi)))+geom_density(position="identity", alpha=0.5, adjust = 1/2) + theme_bw(base_size = 20) + xlab("Vdem_Media") +
  scale_fill_discrete(name = "Regional Election", labels = c("No", "Yes")) +  theme(legend.position = "bottom") + xlab("Media Freedom")
ggsave(filename = "plot/msfs_bi_regional.png",  width = 6, height = 5, dpi = 300, units = "in", device='png') # save last plot

## Figure A3

# density plot of local government index - media freedom 

cor.test(df$MSF, df$v2xel_locelec, 
         method = "pearson")

df %>%
  ggplot(aes(x = MSF, y = v2xel_locelec))+ geom_point(alpha=0.4) + theme_bw(base_size = 20) + #geom_density_2d() +
  xlab('Media Freedom') + ylab('Local Government Index') + annotate("text", x = 0.7, y = 1, label = "Pearson's r = 0.344", col='blue')
ggsave(filename = "plot/msfs_local.png",  width = 6, height = 5, dpi = 500, units = "in", device='png') # save last plot

# density plot of regional government index - media freedom 

cor.test(df$MSF, df$v2xel_regelec, 
         method = "pearson")

df %>%
  ggplot(aes(x = MSF, y = v2xel_regelec))+ geom_point(alpha=0.4) + theme_bw(base_size = 20) + #geom_density_2d() +
  xlab('Media Freedom') + ylab('Regional Government Index') + annotate("text", x = 0.7, y = 1, label = "Pearson's r = 0.334", col='blue')
ggsave(filename = "plot/msfs_regional.png",  width = 6, height = 5, dpi = 500, units = "in", device='png') # save last plot


## Figure A4

# relationship between national elections and local elections

cor.test(df$lag_natelec, df$lag_ellocelc_bi, 
         method = "pearson")

df %>%
  select(lag_natelec, lag_ellocelc_bi) %>%
  drop_na() %>%
  ggplot(aes(x=factor(lag_natelec), y= factor(lag_ellocelc_bi))) + geom_jitter(alpha = .6)+
  theme_bw(base_size = 20) + 
  xlab('National Elections') + ylab('Local Elections')
ggsave(filename = "plot/nat_loc.png",  width = 5, height = 5, dpi = 300, units = "in", device='png') # save last plot

# relationship between national elections and local elections

cor.test(df$lag_natelec, df$v2elsrgel_bi, 
         method = "pearson")

df %>%
  select(lag_natelec, v2elsrgel_bi) %>%
  drop_na() %>%
  ggplot(aes(x=factor(lag_natelec), y= factor(v2elsrgel_bi))) + geom_jitter(alpha = .6)+
  theme_bw(base_size = 20) + 
  xlab('National Elections') + ylab('Regional Elections')
ggsave(filename = "plot/nat_reg.png",  width = 5, height = 5, dpi = 300, units = "in", device='png') # save last plot


#========================#
# Descriptive Statistics #
#========================#

data <- df %>%
  select(noncivilfail, lag_ellocelc_bi, lag_elsrgel_bi, lag_xellocelec, lag_xelregelec,
         lag_msfs, lag_personalism, lag_natelec, lag_gdp, lag_growth, lag_logoil, lag_war, lag_interwar)

var_name <- c("Regime Breakdown", "Local Election", "Regional Election", "Local Quality", "Regional Quality",
              "Media Freedom", "Personalism", "National Election", "ln(GDP)", "Growth", "ln(oil)", "Civil War", "Interstate War")

summary <- tibble(var_name, sapply(data, mean, na.rm=TRUE),
                  sapply(data, sd, na.rm=TRUE),
                  sapply(data, min, na.rm=TRUE),
                  sapply(data, max, na.rm=TRUE)) %>%
  rename(Mean = `sapply(data, mean, na.rm = TRUE)`,
         SD = `sapply(data, sd, na.rm = TRUE)`,
         Min = `sapply(data, min, na.rm = TRUE)`,
         Max = `sapply(data, max, na.rm = TRUE)`)

xtable(summary)

#=================#
# List of regimes #
#=================#

df_group <- df %>%
  select(country, caseid, year) %>%
  group_by(country, caseid) %>%
  slice(c(1,n())) %>%
  summarize(Year = paste0(year, collapse = "-")) %>%
  select(-caseid)

write_excel_csv(df_group, 'data/clean/list_regime.csv')


#==================================================================#
# Modified binary indicators using 10% quantile of quality indices #
#==================================================================#


q10_xellocelec = quantile(df$lag_xellocelec, probs = c(0.1), na.rm = T)
q10_xelregelec = quantile(df$lag_xelregelec, probs = c(0.1), na.rm = T)


df <- read_csv('data/clean/df.csv') %>%
  mutate(lag_ellocelc_q10 = ifelse(lag_ellocelc_bi == 1 & lag_xellocelec <= q10_xellocelec, 0, lag_ellocelc_bi),
         lag_elsrgel_q10 = ifelse(lag_elsrgel_bi == 1 & lag_xelregelec <= q10_xelregelec, 0, lag_elsrgel_bi))


# lag_ellocelc_q10

fe_loc_base <- felm(gwf_case_fail ~ lag_ellocelc_q10 * lag_msfs + 
                      poly(gwf_case_duration,3)+ lag_gdp + lag_growth+
                      lag_logoil + lag_war + lag_natelec + lag_personalism| caseid + year | 0 | caseid,
                    data=df)
summary(fe_loc_base)

# lag_elsrgel_q10

fe_reg_base <- felm(gwf_case_fail ~ lag_elsrgel_q10 * lag_msfs + 
                      poly(gwf_case_duration,3)+ lag_gdp + lag_growth+
                      lag_logoil + lag_war + lag_natelec + lag_personalism| caseid + year | 0 | caseid,
                    data=df)

summary(fe_reg_base)

# report: Table A1
stargazer(fe_loc_base, fe_reg_base, omit.stat=c("rsq", "adj.rsq", "ser"),
          omit = c("lag_gdp", "lag_growth", "lag_logoil", "lag_war", "lag_interwar", 
                   "lag_natelec", "lag_personalism", "gwf_case_duration"),
          covariate.labels=c("Local Election","Regional Election",
                             "Media Freedom", "Local Election * Media Freedm", "Regional Election * Media Freedom"),
          dep.var.labels.include = FALSE, model.numbers = FALSE,
          df=F, column.labels = c("Local","Regional"), 
          dep.var.caption  = "DV: Regime Breakdown",
          add.lines = list(c("Regime-Fixed", "Y", "Y"),
                           c("Year-Fixed", "Y", "Y"),
                           c("Confounders", "Y", "Y"),
                           c("Poly(Duration,3)", "Y", "Y")),
          notes = "SE clustered at the regime in parentheses.",
          title = "Linear Probability Model")

#==================#
# All case failure #
#==================#

# lag_ellocelc_bi

fe_loc_base <- felm(gwf_case_fail ~ lag_ellocelc_bi * lag_msfs + 
                      poly(gwf_case_duration,3)+ lag_gdp + lag_growth+
                      lag_logoil + lag_war + lag_natelec + lag_personalism| caseid + year | 0 | caseid,
                    data=df)


# lag_elsrgel_bi

fe_reg_base <- felm(gwf_case_fail ~ lag_elsrgel_bi * lag_msfs + 
                      poly(gwf_case_duration,3)+ lag_gdp + lag_growth+
                      lag_logoil + lag_war + lag_natelec + lag_personalism| caseid + year | 0 | caseid,
                    data=df)


# lag_xellocelec

fe_locq_base <- felm(gwf_case_fail ~ lag_xellocelec * lag_msfs + 
                       poly(gwf_case_duration,3)+ lag_gdp + lag_growth+
                       lag_logoil + lag_war + lag_natelec + lag_personalism| caseid + year | 0 | caseid,
                     data=df)

# lag_xelregelec

fe_regq_base <- felm(gwf_case_fail ~ lag_xelregelec * lag_msfs + 
                       poly(gwf_case_duration,3)+ lag_gdp + lag_growth+
                       lag_logoil + lag_war + lag_natelec + lag_personalism| caseid + year | 0 | caseid,
                     data=df)

## report: Table A2
stargazer(fe_loc_base, fe_reg_base, fe_locq_base, fe_regq_base, omit.stat=c("rsq", "adj.rsq", "ser"),
          omit = c("lag_gdp", "lag_growth", "lag_logoil", "lag_war", "lag_interwar", 
                   "lag_natelec", "lag_personalism", "gwf_case_duration"),
          covariate.labels=c("Local Election","Regional Election","Local Quality", "Regional Quality",
                             "Media Freedom", "Local Election * Media Freedm", "Regional Election * Media Freedom",
                             "Local Quality * Media Freedom", "Regional Quality * Media Freedom"),
          dep.var.labels.include = FALSE, model.numbers = FALSE,
          df=F, column.labels = c("Local","Regional","Loc Quality", "Reg Quality"), 
          dep.var.caption  = "DV: Regime Breakdown",
          add.lines = list(c("Regime-Fixed", "Y", "Y","Y", "Y"),
                           c("Year-Fixed", "Y", "Y","Y", "Y"),
                           c("Confounders", "Y", "Y","Y", "Y"),
                           c("Poly(Duration,3)", "Y", "Y","Y", "Y")),
          notes = "SE clustered at the regime in parentheses.",
          title = "Linear Probability Model")

#================#
# No interaction #
#================#

# lag_ellocelc_bi

fe_loc_base <- felm(noncivilfail ~ lag_ellocelc_bi + lag_msfs + 
                      poly(gwf_case_duration,3)+ lag_gdp + lag_growth+
                      lag_logoil + lag_war + lag_interwar+lag_natelec + lag_personalism| caseid + year | 0 | caseid,
                    data=df)

# lag_elsrgel_bi

fe_reg_base <- felm(noncivilfail ~ lag_elsrgel_bi + lag_msfs + 
                      poly(gwf_case_duration,3)+ lag_gdp + lag_growth+
                      lag_logoil + lag_war + lag_interwar+lag_natelec + lag_personalism| caseid + year | 0 | caseid,
                    data=df)


# lag_xellocelec

fe_locq_base <- felm(noncivilfail ~ lag_xellocelec + lag_msfs + 
                       poly(gwf_case_duration,3)+ lag_gdp + lag_growth+
                       lag_logoil + lag_war + lag_interwar+lag_natelec + lag_personalism| caseid + year | 0 | caseid,
                     data=df)

# lag_xelregelec

fe_regq_base <- felm(noncivilfail ~ lag_xelregelec + lag_msfs + 
                       poly(gwf_case_duration,3)+ lag_gdp + lag_growth+
                       lag_logoil + lag_war + lag_interwar+lag_natelec + lag_personalism| caseid + year | 0 | caseid,
                     data=df)

## report: Table A3
stargazer(fe_loc_base, fe_reg_base, fe_locq_base, fe_regq_base, omit.stat=c("rsq", "adj.rsq", "ser"),
          omit = c("lag_gdp", "lag_growth", "lag_logoil", "lag_war", "lag_interwar", 
                   "lag_natelec", "lag_personalism", "gwf_case_duration"),
          covariate.labels=c("Local Election","Regional Election","Local Quality", "Regional Quality",
                             "Media Freedom"),
          dep.var.labels.include = FALSE, model.numbers = FALSE,
          df=F, column.labels = c("Local","Regional","Loc Quality", "Reg Quality"), 
          dep.var.caption  = "DV: Regime Breakdown",
          add.lines = list(c("Regime-Fixed", "Y", "Y","Y", "Y"),
                           c("Year-Fixed", "Y", "Y","Y", "Y"),
                           c("Confounders", "Y", "Y","Y", "Y"),
                           c("Poly(Duration,3)", "Y", "Y","Y", "Y")),
          notes = "SE clustered at the regime in parentheses.",
          title = "Linear Probability Model")

#========================================#
# National election interaction included #
#========================================#

# lag_ellocelc_bi

fe_loc_base <- felm(noncivilfail ~ lag_ellocelc_bi * lag_msfs + lag_natelec * lag_msfs + 
                      poly(gwf_case_duration,3)+ lag_gdp + lag_growth+
                      lag_logoil + lag_war + lag_interwar + lag_personalism| caseid + year | 0 | caseid,
                    data=df)

# lag_elsrgel_bi

fe_reg_base <- felm(noncivilfail ~ lag_elsrgel_bi * lag_msfs + lag_natelec * lag_msfs + 
                      poly(gwf_case_duration,3)+ lag_gdp + lag_growth+
                      lag_logoil + lag_war + lag_interwar + lag_personalism| caseid + year | 0 | caseid,
                    data=df)


# lag_xellocelec

fe_locq_base <- felm(noncivilfail ~ lag_xellocelec * lag_msfs + lag_natelec * lag_msfs + 
                       poly(gwf_case_duration,3)+ lag_gdp + lag_growth+
                       lag_logoil + lag_war + lag_interwar + lag_personalism| caseid + year | 0 | caseid,
                     data=df)

# lag_xelregelec

fe_regq_base <- felm(noncivilfail ~ lag_xelregelec * lag_msfs + lag_natelec * lag_msfs + 
                       poly(gwf_case_duration,3)+ lag_gdp + lag_growth+
                       lag_logoil + lag_war + lag_interwar + lag_personalism| caseid + year | 0 | caseid,
                     data=df)

## report: Table A4

stargazer(fe_loc_base, fe_reg_base, fe_locq_base, fe_regq_base, omit.stat=c("rsq", "adj.rsq", "ser"),
          omit = c("lag_gdp", "lag_growth", "lag_logoil", "lag_war", "lag_interwar", 
                  "lag_personalism", "gwf_case_duration"),
          covariate.labels=c("Local Election","Regional Election","Local Quality", "Regional Quality",
                             "Media Freedom", "National Election", 
                             "Local Election * Media Freedm", "Regional Election * Media Freedom",
                             "Local Quality * Media Freedom", "Regional Quality * Media Freedom",
                             "National Election * Media Freedom"),
          dep.var.labels.include = FALSE, model.numbers = FALSE,
          df=F, column.labels = c("Local","Regional","Loc Quality", "Reg Quality"), 
          dep.var.caption  = "DV: Regime Breakdown",
          add.lines = list(c("Regime-Fixed", "Y", "Y","Y", "Y"),
                           c("Year-Fixed", "Y", "Y","Y", "Y"),
                           c("Confounders", "Y", "Y","Y", "Y"),
                           c("Poly(Duration,3)", "Y", "Y","Y", "Y")),
          notes = "SE clustered at the regime in parentheses.",
          title = "Linear Probability Model")


#====================================#
# Fariss's human rights as moderator #
#====================================#

# lag_ellocelc_bi

fe_loc_base <- felm(noncivilfail ~ lag_ellocelc_bi * lag_fariss + 
                      poly(gwf_case_duration,3)+ lag_gdp + lag_growth+
                      lag_logoil + lag_war + lag_interwar+lag_natelec + lag_personalism| caseid + year | 0 | caseid,
                    data=df)

# lag_elsrgel_bi

fe_reg_base <- felm(noncivilfail ~ lag_elsrgel_bi * lag_fariss + 
                      poly(gwf_case_duration,3)+ lag_gdp + lag_growth+
                      lag_logoil + lag_war + lag_interwar+lag_natelec + lag_personalism| caseid + year | 0 | caseid,
                    data=df)


# lag_xellocelec

fe_locq_base <- felm(noncivilfail ~ lag_xellocelec * lag_fariss + 
                       poly(gwf_case_duration,3)+ lag_gdp + lag_growth+
                       lag_logoil + lag_war + lag_interwar+lag_natelec + lag_personalism| caseid + year | 0 | caseid,
                     data=df)

# lag_xelregelec

fe_regq_base <- felm(noncivilfail ~ lag_xelregelec * lag_fariss + 
                       poly(gwf_case_duration,3)+ lag_gdp + lag_growth+
                       lag_logoil + lag_war + lag_interwar+lag_natelec + lag_personalism| caseid + year | 0 | caseid,
                     data=df)

## report: Table A5
stargazer(fe_loc_base, fe_reg_base, fe_locq_base, fe_regq_base, omit.stat=c("rsq", "adj.rsq", "ser"),
          omit = c("lag_gdp", "lag_growth", "lag_logoil", "lag_war", "lag_interwar", 
                   "lag_natelec", "lag_personalism", "gwf_case_duration"),
          covariate.labels=c("Local Election","Regional Election","Local Quality", "Regional Quality",
                             "Human Rights", "Local Election * Human Rights", "Regional Election * Human Rights",
                             "Local Quality * Human Rights", "Regional Quality * Human Rights"),
          dep.var.labels.include = FALSE, model.numbers = FALSE,
          df=F, column.labels = c("Local","Regional","Loc Quality", "Reg Quality"), 
          dep.var.caption  = "DV: Regime Breakdown",
          add.lines = list(c("Regime-Fixed", "Y", "Y","Y", "Y"),
                           c("Year-Fixed", "Y", "Y","Y", "Y"),
                           c("Confounders", "Y", "Y","Y", "Y"),
                           c("Poly(Duration,3)", "Y", "Y","Y", "Y")),
          notes = "SE clustered at the regime in parentheses.",
          title = "Linear Probability Model")

#======================================#
# V-Dem Liberal democracy as moderator #
#======================================#

# lag_ellocelc_bi

fe_loc_base <- felm(noncivilfail ~ lag_ellocelc_bi * lag_libdem + 
                      poly(gwf_case_duration,3)+ lag_gdp + lag_growth+
                      lag_logoil + lag_war + lag_interwar+lag_natelec + lag_personalism| caseid + year | 0 | caseid,
                    data=df)

# lag_elsrgel_bi

fe_reg_base <- felm(noncivilfail ~ lag_elsrgel_bi * lag_libdem + 
                      poly(gwf_case_duration,3)+ lag_gdp + lag_growth+
                      lag_logoil + lag_war + lag_interwar+lag_natelec + lag_personalism| caseid + year | 0 | caseid,
                    data=df)


# lag_xellocelec

fe_locq_base <- felm(noncivilfail ~ lag_xellocelec * lag_libdem + 
                       poly(gwf_case_duration,3)+ lag_gdp + lag_growth+
                       lag_logoil + lag_war + lag_interwar+lag_natelec + lag_personalism| caseid + year | 0 | caseid,
                     data=df)

# lag_xelregelec

fe_regq_base <- felm(noncivilfail ~ lag_xelregelec * lag_libdem + 
                       poly(gwf_case_duration,3)+ lag_gdp + lag_growth+
                       lag_logoil + lag_war + lag_interwar+lag_natelec + lag_personalism| caseid + year | 0 | caseid,
                     data=df)

## report: Table A6
stargazer(fe_loc_base, fe_reg_base, fe_locq_base, fe_regq_base, omit.stat=c("rsq", "adj.rsq", "ser"),
          omit = c("lag_gdp", "lag_growth", "lag_logoil", "lag_war", "lag_interwar", 
                   "lag_natelec", "lag_personalism", "gwf_case_duration"),
          covariate.labels=c("Local Election","Regional Election","Local Quality", "Regional Quality",
                             "Liberalization", "Local Election * Liberalization", "Regional Election * Liberalization",
                             "Local Quality * Liberalization", "Regional Quality * Liberalization"),
          dep.var.labels.include = FALSE, model.numbers = FALSE,
          df=F, column.labels = c("Local","Regional","Loc Quality", "Reg Quality"), 
          dep.var.caption  = "DV: Regime Breakdown",
          add.lines = list(c("Regime-Fixed", "Y", "Y","Y", "Y"),
                           c("Year-Fixed", "Y", "Y","Y", "Y"),
                           c("Confounders", "Y", "Y","Y", "Y"),
                           c("Poly(Duration,3)", "Y", "Y","Y", "Y")),
          notes = "SE clustered at the regime in parentheses.",
          title = "Linear Probability Model")



#=====================#
# V-Dem polity2 as DV #
#=====================#

# lag_ellocelc_bi

fe_loc_base <- felm(e_polity2 ~ lag_ellocelc_bi * lag_msfs + 
                      poly(gwf_case_duration,3)+ lag_gdp + lag_growth+
                      lag_logoil + lag_war + lag_interwar+lag_natelec + lag_personalism| caseid + year | 0 | caseid,
                    data=df)

# lag_elsrgel_bi

fe_reg_base <- felm(e_polity2 ~ lag_elsrgel_bi * lag_msfs + 
                      poly(gwf_case_duration,3)+ lag_gdp + lag_growth+
                      lag_logoil + lag_war + lag_interwar+lag_natelec + lag_personalism| caseid + year | 0 | caseid,
                    data=df)


# lag_xellocelec

fe_locq_base <- felm(e_polity2 ~ lag_xellocelec * lag_msfs + 
                       poly(gwf_case_duration,3)+ lag_gdp + lag_growth+
                       lag_logoil + lag_war + lag_interwar+lag_natelec + lag_personalism| caseid + year | 0 | caseid,
                     data=df)

# lag_xelregelec

fe_regq_base <- felm(e_polity2 ~ lag_xelregelec * lag_msfs + 
                       poly(gwf_case_duration,3)+ lag_gdp + lag_growth+
                       lag_logoil + lag_war + lag_interwar+lag_natelec + lag_personalism| caseid + year | 0 | caseid,
                     data=df)

## report: Table A7
stargazer(fe_loc_base, fe_reg_base, fe_locq_base, fe_regq_base, omit.stat=c("rsq", "adj.rsq", "ser"),
          omit = c("lag_gdp", "lag_growth", "lag_logoil", "lag_war", "lag_interwar", 
                   "lag_natelec", "lag_personalism", "gwf_case_duration"),
          covariate.labels=c("Local Election","Regional Election","Local Quality", "Regional Quality",
                             "Media Freedom", "Local Election * Media Freedom", "Regional Election * Media Freedom",
                             "Local Quality * Media Freedom", "Regional Quality * Media Freedom"),
          dep.var.labels.include = FALSE, model.numbers = FALSE,
          df=F, column.labels = c("Local","Regional","Loc Quality", "Reg Quality"), 
          dep.var.caption  = "DV: Polity2",
          add.lines = list(c("Regime-Fixed", "Y", "Y","Y", "Y"),
                           c("Year-Fixed", "Y", "Y","Y", "Y"),
                           c("Confounders", "Y", "Y","Y", "Y"),
                           c("Poly(Duration,3)", "Y", "Y","Y", "Y")),
          notes = "SE clustered at the regime in parentheses.",
          title = "Linear Regression Model")



df %>%
  filter(country == "Brazil") %>%
  ggplot(aes(x= year, y = e_polity2)) + geom_point() + geom_line()


#===========================================================#
# Analaysis with country-year with national-level elections #
#===========================================================#


df <- read_csv('data/clean/df.csv') %>%
  filter(lag_natelec == 1)

# lag_ellocelc_bi

fe_loc_base <- felm(noncivilfail ~ lag_ellocelc_bi * lag_msfs + 
                      poly(gwf_case_duration,3)+ lag_gdp + lag_growth+
                      lag_logoil + lag_war  + lag_personalism| caseid + year | 0 | caseid,
                    data=df)

# lag_elsrgel_bi

fe_reg_base <- felm(noncivilfail ~ lag_elsrgel_bi * lag_msfs + 
                      poly(gwf_case_duration,3)+ lag_gdp + lag_growth+
                      lag_logoil + lag_war  + lag_personalism| caseid + year | 0 | caseid,
                    data=df)


# lag_xellocelec

fe_locq_base <- felm(noncivilfail ~ lag_xellocelec * lag_msfs + 
                       poly(gwf_case_duration,3)+ lag_gdp + lag_growth+
                       lag_logoil + lag_war  + lag_personalism| caseid + year | 0 | caseid,
                     data=df)

# lag_xelregelec

fe_regq_base <- felm(noncivilfail ~ lag_xelregelec * lag_msfs + 
                       poly(gwf_case_duration,3)+ lag_gdp + lag_growth+
                       lag_logoil + lag_war  + lag_personalism| caseid + year | 0 | caseid,
                     data=df)

## report: Table A9
stargazer(fe_loc_base, fe_reg_base, fe_locq_base, fe_regq_base, omit.stat=c("rsq", "adj.rsq", "ser"),
          omit = c("lag_gdp", "lag_growth", "lag_logoil", "lag_war", "lag_interwar", 
                   "lag_natelec", "lag_personalism", "gwf_case_duration"),
          covariate.labels=c("Local Election","Regional Election","Local Quality", "Regional Quality",
                             "Media Freedom", "Local Election * Media Freedom", "Regional Election * Media Freedom",
                             "Local Quality * Media Freedom", "Regional Quality * Media Freedom"),
          dep.var.labels.include = FALSE, model.numbers = FALSE,
          df=F, column.labels = c("Local","Regional","Loc Quality", "Reg Quality"), 
          dep.var.caption  = "DV: Regime Breakdown",
          add.lines = list(c("Regime-Fixed", "Y", "Y","Y", "Y"),
                           c("Year-Fixed", "Y", "Y","Y", "Y"),
                           c("Confounders", "Y", "Y","Y", "Y"),
                           c("Poly(Duration,3)", "Y", "Y","Y", "Y")),
          notes = "SE clustered at the regime in parentheses.",
          title = "Linear Probability Model")



#====================================#
# Multiple Imputation using Amelia 2 #
#====================================#

setwd('C:/Users/battl/OneDrive/JUN/UC MERCED/Independent Research/local election and media freedom')
df_imp <- read_csv('data/clean/df.csv') %>%
  select(cowcode, year, lag_msfs, gwf_case_duration, lag_gdp, lag_growth, 
         lag_logoil, lag_personalism)

df_imp <- as.data.frame(df_imp)

set.seed(123)
a.out <- amelia(df_imp, m = 5, ts = "year", cs = "cowcode")

setwd('C:/Users/battl/OneDrive/JUN/UC MERCED/Independent Research/local election and media freedom/data/imputation')
write.amelia(obj=a.out, file.stem = "outdata")

df_mainv <- df %>%
  select(cowcode, year, noncivilfail, lag_ellocelc_bi, lag_elsrgel_bi, 
         lag_xellocelec, lag_xelregelec, lag_war, lag_interwar, lag_natelec, caseid)

df_imp <- read_csv('outdata2.csv') %>%
  left_join(df_mainv) %>%
  select(-1)


write_excel_csv(df_imp, "C:/Users/battl/OneDrive/JUN/UC MERCED/Independent Research/local election and media freedom/data/clean/imputation.csv")


# lag_ellocelc_bi

fe_loc_base <- felm(noncivilfail ~ lag_ellocelc_bi * lag_msfs + 
                      poly(gwf_case_duration,3)+ lag_gdp + lag_growth+
                      lag_logoil + lag_war + lag_interwar+lag_natelec + lag_personalism| caseid + year | 0 | caseid,
                    data=df_imp)

# lag_elsrgel_bi

fe_reg_base <- felm(noncivilfail ~ lag_elsrgel_bi * lag_msfs + 
                      poly(gwf_case_duration,3)+ lag_gdp + lag_growth+
                      lag_logoil + lag_war + lag_interwar+lag_natelec + lag_personalism| caseid + year | 0 | caseid,
                    data=df_imp)


# lag_xellocelec

fe_locq_base <- felm(noncivilfail ~ lag_xellocelec * lag_msfs + 
                       poly(gwf_case_duration,3)+ lag_gdp + lag_growth+
                       lag_logoil + lag_war + lag_interwar+lag_natelec + lag_personalism| caseid + year | 0 | caseid,
                     data=df_imp)

# lag_xelregelec

fe_regq_base <- felm(noncivilfail ~ lag_xelregelec * lag_msfs + 
                       poly(gwf_case_duration,3)+ lag_gdp + lag_growth+
                       lag_logoil + lag_war + lag_interwar+lag_natelec + lag_personalism| caseid + year | 0 | caseid,
                     data=df_imp)

## report: Table A10
stargazer(fe_loc_base, fe_reg_base, fe_locq_base, fe_regq_base, omit.stat=c("rsq", "adj.rsq", "ser"),
          omit = c("lag_gdp", "lag_growth", "lag_logoil", "lag_war", "lag_interwar", 
                   "lag_natelec", "lag_personalism", "gwf_case_duration"),
          covariate.labels=c("Local Election","Regional Election","Local Quality", "Regional Quality",
                             "Media Freedom", "Local Election * Media Freedom", "Regional Election * Media Freedom",
                             "Local Quality * Media Freedom", "Regional Quality * Media Freedom"),
          dep.var.labels.include = FALSE, model.numbers = FALSE,
          df=F, column.labels = c("Local","Regional","Loc Quality", "Reg Quality"), 
          dep.var.caption  = "DV: Regime Breakdown",
          add.lines = list(c("Regime-Fixed", "Y", "Y","Y", "Y"),
                           c("Year-Fixed", "Y", "Y","Y", "Y"),
                           c("Confounders", "Y", "Y","Y", "Y"),
                           c("Poly(Duration,3)", "Y", "Y","Y", "Y")),
          notes = "SE clustered at the regime in parentheses.",
          title = "Linear Probability Model")


#===================#
# CID fixed effects #
#===================#

# lag_ellocelc_bi

fe_loc_base <- felm(noncivilfail ~ lag_ellocelc_bi * lag_msfs + 
                      poly(gwf_case_duration,3)+ lag_gdp + lag_growth+
                      lag_logoil + lag_war + lag_interwar+lag_natelec + gwf_military + gwf_personal + gwf_party| cid + year | 0 | cid,
                    data=df)

# lag_elsrgel_bi

fe_reg_base <- felm(noncivilfail ~ lag_elsrgel_bi * lag_msfs + 
                      poly(gwf_case_duration,3)+ lag_gdp + lag_growth+
                      lag_logoil + lag_war + lag_interwar+lag_natelec + gwf_military + gwf_personal + gwf_party| cid + year | 0 | cid,
                    data=df)


# lag_xellocelec

fe_locq_base <- felm(noncivilfail ~ lag_xellocelec * lag_msfs + 
                       poly(gwf_case_duration,3)+ lag_gdp + lag_growth+
                       lag_logoil + lag_war + lag_interwar+lag_natelec + gwf_military + gwf_personal + gwf_party| cid + year | 0 | cid,
                     data=df)

# lag_xelregelec

fe_regq_base <- felm(noncivilfail ~ lag_xelregelec * lag_msfs + 
                       poly(gwf_case_duration,3)+ lag_gdp + lag_growth+
                       lag_logoil + lag_war + lag_interwar+lag_natelec + gwf_military + gwf_personal + gwf_party| cid + year | 0 | cid,
                     data=df)

# lag_natelec

fe_nat <- felm(noncivilfail ~ lag_natelec * lag_msfs + 
                 poly(gwf_case_duration,3)+ lag_gdp + lag_growth+
                 lag_logoil + lag_war + lag_interwar+ gwf_military + gwf_personal + gwf_party| cid + year | 0 | cid,
               data=df)


liblocbi <- felm(noncivilfail ~ lag_ellocelc_bi * lag_libdem + 
                   poly(gwf_case_duration,3)+ lag_gdp + lag_growth+
                   lag_logoil + lag_war + lag_interwar+lag_natelec + gwf_military + gwf_personal + gwf_party| cid + year | 0 | cid,
                 data=df)

libregbi <- felm(noncivilfail ~ lag_elsrgel_bi * lag_libdem + 
                   poly(gwf_case_duration,3)+ lag_gdp + lag_growth+
                   lag_logoil + lag_war + lag_interwar+lag_natelec + gwf_military + gwf_personal + gwf_party| cid + year | 0 | cid,
                 data=df)

libloq <- felm(noncivilfail ~ lag_xellocelec * lag_libdem + 
                 poly(gwf_case_duration,3)+ lag_gdp + lag_growth+
                 lag_logoil + lag_war + lag_interwar+lag_natelec + gwf_military + gwf_personal + gwf_party| cid + year | 0 | cid,
               data=df)

libreq <- felm(noncivilfail ~ lag_xelregelec * lag_libdem + 
                 poly(gwf_case_duration,3)+ lag_gdp + lag_growth+
                 lag_logoil + lag_war + lag_interwar+lag_natelec + gwf_military + gwf_personal + gwf_party| cid + year | 0 | cid,
               data=df)


## report: Table A11
stargazer(fe_loc_base, fe_reg_base, fe_locq_base, fe_regq_base, fe_nat,
          liblocbi, libregbi, libloq, libreq,
          omit.stat=c("rsq", "adj.rsq", "ser"),
          omit = c("lag_gdp", "lag_growth", "lag_logoil", "lag_war", "lag_interwar", 
                   "gwf_case_duration", "gwf_military", "gwf_personal", "gwf_party"),
          covariate.labels=c("Local Election","Regional Election","Local Quality", "Regional Quality",
                             "Media Freedom","Liberalization","National Election", "Local Election * Media Freedom", "Regional Election * Media Freedom",
                             "Local Quality * Media Freedom", "Regional Quality * Media Freedom", "National Election * Media Freedom",
                             "Local Election * Liberalization", "Regional Election * Liberalization",
                             "Local Quality * Liberalization", "Regional Quality * Liberalization"),
          dep.var.labels.include = FALSE, model.numbers = FALSE,
          df=F, column.labels = c("Local","Regional","Loc Quality", "Reg Quality", "National Elec"), 
          dep.var.caption  = "DV: Regime Breakdown",
          add.lines = list(c("Country-Fixed", "Y", "Y","Y", "Y","Y"),
                           c("Year-Fixed", "Y", "Y","Y", "Y", "Y"),
                           c("Confounders", "Y", "Y","Y", "Y", "Y"),
                           c("Poly(Duration,3)", "Y", "Y","Y", "Y", "Y")),
          notes = "SE clustered at the country in parentheses.",
          title = "Linear Probability Model")

